
R version 3.2.5 (2016-04-14) -- "Very, Very Secure Dishes"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

Microsoft R Open 3.2.5
Default CRAN mirror snapshot taken on 2016-05-01
The enhanced R distribution from Microsoft
Visit https://mran.microsoft.com/ for information
about additional features.

Multithreaded BLAS/LAPACK libraries detected. Using 3 cores for math algorithms.
[Previously saved workspace restored]

> #################
> # Lebovic and Voeten Replication
> # with Causal Endogeneity Model
> #################
> 
> rm(list=ls())
> set.seed(123456)
> library(foreign)
> library(lme4)
Loading required package: Matrix
> library(mvtnorm)
> dat<-read.dta("jprworkdatanewMOD.dta")
> 
> # the Lebovic and Voeten replication model
> col1<-lmer(BIPOP~l_BIPOP+PUBRES+d_HRIGHTS+l_HRIGHTS+d_CIVIL+l_CIVIL+l_GDPPOP+l_LNPOP+l_USAGREE+WAR+CAPAB+linear+quad+(1 | CCODE), data=dat, REML=F)
Warning message:
Some predictor variables are on very different scales: consider rescaling 
> options(scipen=100)
> summary(col1)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: BIPOP ~ l_BIPOP + PUBRES + d_HRIGHTS + l_HRIGHTS + d_CIVIL +  
    l_CIVIL + l_GDPPOP + l_LNPOP + l_USAGREE + WAR + CAPAB +  
    linear + quad + (1 | CCODE)
   Data: dat

     AIC      BIC   logLik deviance df.resid 
  3613.4   3705.4  -1790.7   3581.4     2308 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.1913 -0.5088 -0.0021  0.5082  6.6941 

Random effects:
 Groups   Name        Variance Std.Dev.
 CCODE    (Intercept) 0.1484   0.3852  
 Residual             0.2407   0.4906  
Number of obs: 2324, groups:  CCODE, 118

Fixed effects:
               Estimate  Std. Error t value
(Intercept)  4.90283014  0.55519972    8.83
l_BIPOP      0.55343082  0.01683916   32.87
PUBRES      -0.02991082  0.05902802   -0.51
d_HRIGHTS   -0.03700406  0.01684928   -2.20
l_HRIGHTS   -0.01752621  0.01817769   -0.96
d_CIVIL     -0.01548797  0.01934723   -0.80
l_CIVIL     -0.02047250  0.01264284   -1.62
l_GDPPOP    -0.10718221  0.02445462   -4.38
l_LNPOP     -0.16223521  0.03035736   -5.34
l_USAGREE   -0.12335610  0.14296053   -0.86
WAR         -0.08765625  0.03438182   -2.55
CAPAB       -0.75853549  0.60532203   -1.25
linear      -0.00491270  0.00109635   -4.48
quad        -0.00004620  0.00008544   -0.54

Correlation of Fixed Effects:
          (Intr) l_BIPO PUBRES d_HRIG l_HRIG d_CIVI l_CIVI l_GDPP l_LNPO l_USAG
l_BIPOP   -0.319                                                               
PUBRES    -0.054  0.045                                                        
d_HRIGHTS -0.014  0.065 -0.119                                                 
l_HRIGHTS  0.003  0.100 -0.218  0.556                                          
d_CIVIL   -0.068  0.049 -0.049 -0.004 -0.115                                   
l_CIVIL   -0.184  0.090 -0.089 -0.049 -0.158  0.363                            
l_GDPPOP  -0.513  0.145  0.070  0.065  0.101  0.025  0.115                     
l_LNPOP   -0.926  0.181  0.067 -0.061 -0.125  0.030  0.047  0.220              
l_USAGREE -0.068  0.017 -0.061  0.008 -0.010  0.030  0.081 -0.085  0.008       
WAR        0.006  0.061 -0.070 -0.221 -0.326  0.020  0.049  0.060 -0.029  0.014
CAPAB      0.477  0.059 -0.004  0.005  0.010 -0.015 -0.040 -0.137 -0.538  0.018
linear     0.229  0.090 -0.102 -0.020 -0.067  0.103  0.169 -0.221 -0.268  0.573
quad      -0.021  0.103  0.056 -0.004  0.003  0.001  0.009  0.078  0.027 -0.688
          WAR    CAPAB  linear
l_BIPOP                       
PUBRES                        
d_HRIGHTS                     
l_HRIGHTS                     
d_CIVIL                       
l_CIVIL                       
l_GDPPOP                      
l_LNPOP                       
l_USAGREE                     
WAR                           
CAPAB      0.017              
linear     0.038  0.142       
quad       0.043 -0.012 -0.488
fit warnings:
Some predictor variables are on very different scales: consider rescaling
> 
> # the state dependence model
> causend<-lmer(BIPOP~l_BIPOP+PUBRES+RESxBI+d_HRIGHTS+l_HRIGHTS+d_CIVIL+l_CIVIL+l_GDPPOP+l_LNPOP+l_USAGREE+WAR+CAPAB+linear+quad+(1 | CCODE), data=dat, REML=F)
Warning message:
Some predictor variables are on very different scales: consider rescaling 
> options(scipen=100)
> summary(causend)
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: BIPOP ~ l_BIPOP + PUBRES + RESxBI + d_HRIGHTS + l_HRIGHTS + d_CIVIL +  
    l_CIVIL + l_GDPPOP + l_LNPOP + l_USAGREE + WAR + CAPAB +  
    linear + quad + (1 | CCODE)
   Data: dat

     AIC      BIC   logLik deviance df.resid 
  3602.5   3700.2  -1784.2   3568.5     2307 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.1298 -0.5166 -0.0119  0.5102  6.6896 

Random effects:
 Groups   Name        Variance Std.Dev.
 CCODE    (Intercept) 0.1421   0.3770  
 Residual             0.2398   0.4897  
Number of obs: 2324, groups:  CCODE, 118

Fixed effects:
               Estimate  Std. Error t value
(Intercept)  4.97119046  0.54724606   9.084
l_BIPOP      0.53934353  0.01738643  31.021
PUBRES      -0.39157695  0.11655444  -3.360
RESxBI       0.13577504  0.03767856   3.604
d_HRIGHTS   -0.03814668  0.01681472  -2.269
l_HRIGHTS   -0.01908175  0.01812750  -1.053
d_CIVIL     -0.01502838  0.01930372  -0.779
l_CIVIL     -0.01932678  0.01259237  -1.535
l_GDPPOP    -0.11726556  0.02431789  -4.822
l_LNPOP     -0.15912419  0.02982965  -5.334
l_USAGREE   -0.14495237  0.14269209  -1.016
WAR         -0.09736689  0.03440429  -2.830
CAPAB       -0.83851619  0.59495851  -1.409
linear      -0.00491655  0.00109094  -4.507
quad        -0.00003770  0.00008528  -0.442

Correlation of Fixed Effects:
          (Intr) l_BIPO PUBRES RESxBI d_HRIG l_HRIG d_CIVI l_CIVI l_GDPP l_LNPO
l_BIPOP   -0.323                                                               
PUBRES    -0.066  0.252                                                        
RESxBI     0.045 -0.267 -0.863                                                 
d_HRIGHTS -0.014  0.068 -0.043 -0.020                                          
l_HRIGHTS  0.003  0.104 -0.086 -0.027  0.556                                   
d_CIVIL   -0.069  0.045 -0.033  0.010 -0.004 -0.115                            
l_CIVIL   -0.185  0.079 -0.073  0.033 -0.050 -0.158  0.362                     
l_GDPPOP  -0.517  0.170  0.123 -0.103  0.066  0.103  0.025  0.114              
l_LNPOP   -0.923  0.171  0.019  0.017 -0.063 -0.128  0.030  0.048  0.217       
l_USAGREE -0.071  0.029  0.013 -0.051  0.009 -0.008  0.029  0.079 -0.079  0.007
WAR        0.003  0.080  0.035 -0.082 -0.219 -0.323  0.019  0.046  0.069 -0.031
CAPAB      0.473  0.068  0.031 -0.038  0.006  0.013 -0.015 -0.041 -0.131 -0.539
linear     0.222  0.089 -0.042 -0.011 -0.019 -0.066  0.103  0.168 -0.216 -0.264
quad      -0.020  0.091  0.003  0.029 -0.004  0.002  0.001  0.010  0.074  0.027
          l_USAG WAR    CAPAB  linear
l_BIPOP                              
PUBRES                               
RESxBI                               
d_HRIGHTS                            
l_HRIGHTS                            
d_CIVIL                              
l_CIVIL                              
l_GDPPOP                             
l_LNPOP                              
l_USAGREE                            
WAR        0.018                     
CAPAB      0.020  0.020              
linear     0.574  0.039  0.141       
quad      -0.688  0.041 -0.013 -0.488
fit warnings:
Some predictor variables are on very different scales: consider rescaling
> 
> # prepare for marginal effects calculations
> beta.draws<-rmvnorm(1000, mean=causend@beta, sigma=as.matrix(vcov(causend)))
> 
> # instantaneous effects of a CHR Resolution
> y.seq<-seq(from=0, to=7.25, by=0.1)
> me<-matrix(data=NA, ncol=3, nrow=length(y.seq))
> 
> for(i in 1:length(y.seq)){
+   
+   y.test<-y.seq[i]
+   me.draws<-beta.draws[,3]+beta.draws[,4]*y.test
+   me[i,]<-quantile(me.draws, probs=c(0.025, 0.5, 0.975))
+   
+   
+ }
> 
> # plot instantaneous marginal effects
> pdf("instant_mfx_applied_RE_r.pdf")
> par(mar=c(5,5,4,2))
> plot(NA, ylim=c(-0.6, 1), xlim=c(0,7.5), xaxt="n", main=c("Instantaneous Marginal Effect","of UNCHR Resolution on Aggregate Bilateral Aid"), xlab=expression(ln(Aid + 1)[it-1]), ylab=expression(partialdiff*ln(Aid + 1)[it]/partialdiff*Resolution[i(t-1)]))
> axis(side=1, at=0:7)
> hist.back<-hist(dat$l_BIPOP, plot=FALSE)
> lines(hist.back, border=hsv(0, 0, 0.5, alpha=1), col=hsv(0, 0, 0.5, alpha=0.25), freq=FALSE, alpha=1)
> 
> lines(me[,2]~y.seq)
> lines(me[,1]~y.seq, lty=2)
> lines(me[,3]~y.seq, lty=2)
> abline(h=0, lty=3)
> 
> legend("topleft", legend=c("marginal effect", "95% confidence interval"), lty=c(1,2))
> text(x=5.75, y=-0.6, "histogram indicates data density", cex=0.75)
> dev.off()
null device 
          1 
> 
> 
> 
> 
> # long term marginal effects of a CHR resolution
> z<-seq(from=0, to=3.25, by=0.05)
> ss.orig<-c()
> ss.hi<-c()
> ss.lo<-c()
> ss.med<-c()
> for(i in 1:length(z)){
+   
+   x.lo<-0
+   x.hi<-1
+   
+   ss.orig[i]<-mean((z[i]+beta.draws[,3]*x.lo)/(1-beta.draws[,2]-beta.draws[,4]*x.lo))
+   ss.draws<-((z[i]+beta.draws[,3]*x.hi)/(1-beta.draws[,2]-beta.draws[,4]*x.hi))-((z[i]+beta.draws[,3]*x.lo)/(1-beta.draws[,2]-beta.draws[,4]*x.lo))
+   ss.quant<-quantile(ss.draws, probs=c(0.025, 0.5, 0.975))
+   ss.hi[i]<-ss.quant[3]
+   ss.lo[i]<-ss.quant[1]
+   ss.med[i]<-ss.quant[2]
+   
+ }
> 
> summary(ss.orig)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   1.766   3.531   3.531   5.297   7.062 
> 
> pdf("longtermME_applied_r.eps")
> plot(ss.med~ss.orig, type="l", ylim=c(-3, 4), main=c("Long Term Marginal Effect", "of UNCHR Resolution on Bilateral Aid"), xlab="ln(Bilateral Aid + 1) steady state, no resolution", ylab="change in ln(Bilateral Aid + 1) steady state after resolution")
> lines(ss.hi~ss.orig, lty=2)
> lines(ss.lo~ss.orig, lty=2)
> abline(h=0, lty=3)
> legend("bottomright", lty=c(1,2), legend=c("marginal effect", "95% confidence interval"))
> dev.off()
null device 
          1 
> 
> 
> 
> 
> 
> 
> ss.hi<-c()
> ss.lo<-c()
> ss.med<-c()
> ss.draws<-matrix(data=NA, ncol=2, nrow=1000)
> 
> 
> # long term effect of a human rights violation
> x<-c(0,1)
> for(i in 1:2){
+   ss.draws[,i]<-(beta.draws[,6]*4)/(1-beta.draws[,2]-beta.draws[,4]*x[i])
+   ss.quant<-quantile(ss.draws[,i], probs=c(0.025, 0.5, 0.975))
+   ss.hi[i]<-ss.quant[3]
+   ss.lo[i]<-ss.quant[1]
+   ss.med[i]<-ss.quant[2]
+ }
> 
> library(gplots)

Attaching package: 'gplots'

The following object is masked from 'package:stats':

    lowess

> pdf("longtermME_Z_applied_r.eps")
> par(mar=c(4, 5, 4, 2))
> plotCI(x=x, y=ss.med, ui=ss.hi, li=ss.lo, ylim=c(-0.85, 0.25), xlim=c(-0.5, 1.5), xaxt="n", yaxt="n", main=c("Long Term Marginal Effect", "of Physical Integrity Abuse on Bilateral Aid"), ylab="", xlab="")
> axis(1, at=0:1, lab=c("w/o UNCHR Resolution","w/ UNCHR Resolution"))
> axis(2, at=seq(from=-0.8, to =0.2, by=0.2))
> mtext(text="change in ln(Bilateral Aid + 1) steady state,", side=2, line=4)
> mtext(text="abuse from min to max", side=2, line=3)
> abline(h=0, lty=3)
> dev.off()
null device 
          1 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> # plot long term trajectories for the modified LV model
> 
> x<-matrix(dat=c(1,1,0,0,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
> y<-rbind(1,x%*%t(beta.draws)-0.75)
> i<-2
> k<-2
> 
> while(round(y[i,1], 7)!=round(y[i-1,1],7)){
+   y<-rbind(y,NA)
+   i<-i+1
+   k<-k+1
+   for(j in 1:1000){
+     x<-matrix(dat=c(1,y[i-1,j],0,y[i-1,j]*0,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
+     y[i,j]<-x%*%t(t(beta.draws[j,]))-0.75
+   }
+   
+ }
> 
> 
> y<-rbind(y,NA)
> i<-i+1
> for(j in 1:1000){
+   x<-matrix(dat=c(1,y[i-1,j],1,y[i-1,j]*1,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
+   y[i,j]<-x%*%t(t(beta.draws[j,]))-0.75
+ }
> 
> pen<-i+9
> 
> while(i<pen){
+   i<-i+1
+   y<-rbind(y,NA)
+   for(j in 1:1000){
+     x<-matrix(dat=c(1,y[i-1,j],1,y[i-1,j]*1,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
+     y[i,j]<-x%*%t(t(beta.draws[j,]))-0.75
+   }
+ }
> 
> while(round(y[i,1], 7)!=round(y[i-1,1],7)){
+   y<-rbind(y,NA)
+   i<-i+1
+   for(j in 1:1000){
+     x<-matrix(dat=c(1,y[i-1,j],0,y[i-1,j]*0,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
+     y[i,j]<-x%*%t(t(beta.draws[j,]))-0.75
+   }
+   
+ }
> 
> end.point<-(k-10)
> t<-0:(i-end.point)
> plot(y[end.point:i]~t, type="l", xlab="time")
> abline(v=10, lty=2)
> abline(v=20, lty=2)
> 
> traj.quant<-apply(X=y, MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
> 
> plot(traj.quant[2,end.point:i]~t, type="l", xlab="time")
> lines(traj.quant[1,end.point:i]~t, lty=2)
> lines(traj.quant[3,end.point:i]~t, lty=2)
> 
> library(gplots)
> pdf("smallaidtrajectory.eps")
> plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=traj.quant[1,end.point:i], ylab="ln(aggregate bilateral aid PC + 1)", xlab="time", ylim=c(.41, 1.81), main = c("simulated aid trajectory with 90% CI error bars"))
> rect(11,-1,20,4,density=NA, col="grey85")
> plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=traj.quant[1,end.point:i], add=T)
> box()
> abline(v=11, lty=2, lwd=2)
> abline(v=20, lty=2, lwd=2)
> legend("bottomright", legend=c(" median aid prediction", " condemnation period"), pt.cex=c(1,2), pch=c(1,22), pt.bg="gray85")
> dev.off()
pdf 
  2 
> 
> traj.quant[2, (end.point+10):(end.point+20)]
 [1] 1.5184587 1.3349906 1.2079921 1.1295782 1.0753371 1.0372243 1.0119462
 [8] 0.9941234 0.9815145 0.9733852 0.9681789
> # percentage change in aid
> 1 - exp(traj.quant[2, (end.point+20)]) / exp(traj.quant[2, (end.point+10)])
      50% 
0.4232116 
> 
> 
> 
> 
> 
> 
> x<-matrix(dat=c(1,1,0,0,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
> y<-rbind(1,x%*%t(beta.draws)+0.5)
> i<-2
> k<-2
> 
> while(round(y[i,1], 7)!=round(y[i-1,1],7)){
+   y<-rbind(y,NA)
+   i<-i+1
+   k<-k+1
+   for(j in 1:1000){
+     x<-matrix(dat=c(1,y[i-1,j],0,y[i-1,j]*0,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
+     y[i,j]<-x%*%t(t(beta.draws[j,]))+0.5
+   }
+   
+ }
> 
> 
> y<-rbind(y,NA)
> i<-i+1
> for(j in 1:1000){
+   x<-matrix(dat=c(1,y[i-1,j],1,y[i-1,j]*1,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
+   y[i,j]<-x%*%t(t(beta.draws[j,]))+0.5
+ }
> 
> pen<-i+9
> 
> while(i<pen){
+   i<-i+1
+   y<-rbind(y,NA)
+   for(j in 1:1000){
+     x<-matrix(dat=c(1,y[i-1,j],1,y[i-1,j]*1,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
+     y[i,j]<-x%*%t(t(beta.draws[j,]))+0.5
+   }
+ }
> 
> while(round(y[i,1], 7)!=round(y[i-1,1],7)){
+   y<-rbind(y,NA)
+   i<-i+1
+   for(j in 1:1000){
+     x<-matrix(dat=c(1,y[i-1,j],0,y[i-1,j]*0,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
+     y[i,j]<-x%*%t(t(beta.draws[j,]))+0.5
+   }
+   
+ }
> 
> end.point<-(k-10)
> 
> t<-0:(i-end.point)
> plot(y[end.point:i,1]~t, type="l", xlab="time")
> abline(v=10, lty=2)
> abline(v=20, lty=2)
> 
> traj.quant<-apply(X=y, MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
> 
> 
> library(gplots)
> pdf("largeaidtrajectory.eps")
> plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=traj.quant[1,end.point:i], ylab="ln(aggregate bilateral aid PC + 1)", xlab="time", main = c("simulated aid trajectory with 90% CI error bars"), ylim=c(3.99, 5.41))
> rect(11,0,20,6,density=NA, col="grey85")
> plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=traj.quant[1,end.point:i], add=T)
> box()
> abline(v=11, lty=2, lwd=2)
> abline(v=20, lty=2, lwd=2)
> legend("topright", legend=c(" median aid prediction", " condemnation period"), pt.cex=c(1,2), pch=c(1,22), pt.bg="gray85")
> dev.off()
pdf 
  2 
> 
> traj.quant[2, (end.point+10):(end.point+20)]
 [1] 4.231679 4.408053 4.531096 4.613492 4.666653 4.705913 4.732538 4.750149
 [9] 4.761840 4.769255 4.774185
> # percentage change in aid
> exp(traj.quant[2, (end.point+20)]) / exp(traj.quant[2, (end.point+10)]) - 1
      50% 
0.7203111 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> ### plot instantaneous MEs, Neilsen Model 3
> 
> 
> rm(list=ls())
> library(MASS)
> require(foreign)
> 
> dat.neilsen<-read.dta(file="neilsen_3_R_figure_data.dta")
> 
> neilsen.beta<-read.delim(file="xttobit_beta.txt", header=T)
> neilsen.beta<-as.matrix(neilsen.beta[2:31])
> names(neilsen.beta)<-NULL
> 
> neilsen.vcv<-read.delim(file="xttobit_VCV.txt")
> neilsen.vcv<-as.matrix(neilsen.vcv[1:30,2:31])
> 
> beta.draws<-mvrnorm(n=1000, mu=neilsen.beta, Sigma=neilsen.vcv)
> 
> 
> # instantaneous effects of a CHR Resolution
> par(mar=c(5,5,4,2))
> y.seq<-seq(from=0, to=13, by=0.1)
> me<-matrix(data=NA, ncol=3, nrow=length(y.seq))
> 
> for(i in 1:length(y.seq)){
+   
+   y.test<-y.seq[i]
+   me.draws<-beta.draws[,1]+beta.draws[,2]*y.test
+   me[i,]<-quantile(me.draws, probs=c(0.025, 0.5, 0.975))
+   
+   
+ }
> 
> pdf("neilsen_instant.pdf")
> # plot instantaneous marginal effects
> plot(NA, ylim=c(-2, 2.5), xlim=c(0,13), xaxt="n", main=c("Instantaneous Marginal Effect","of UNCHR Resolution on Bilateral Economic Aid"), xlab=expression(ln("Bilateral Economic Aid + 1")[it-1]), ylab=expression(partialdiff*ln("Bilateral Economic Aid + 1")[it]/partialdiff*Resolution[i(t-1)]))
> axis(side=1, at=seq(from=0, to=12, by=2))
> hist.back<-hist(dat.neilsen$l_lneconaidpc, plot=FALSE)
> hist.back$density <- hist.back$density*3
> lines(hist.back, border=hsv(0, 0, 0.5, alpha=1), col=hsv(0, 0, 0.5, alpha=0.25), freq=FALSE, alpha=1)
> 
> lines(me[,2]~y.seq)
> lines(me[,1]~y.seq, lty=2)
> lines(me[,3]~y.seq, lty=2)
> abline(h=0, lty=3)
> 
> legend("topleft", legend=c("marginal effect", "95% confidence interval"), lty=c(1,2))
> text(x=10, y=-2, "histogram indicates 3*(data density)", cex=0.75)
> dev.off()
pdf 
  2 
> 
> 
> ### plot long term trajectories, Neilsen model 3
> 
> x<-matrix(dat=c(0, 2*0, 3.880258, 0, 0*0, 0, 0*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, 2, 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
> y<-rbind(1,x%*%t(beta.draws)+8.8)
> i<-2
> k<-2
> 
> while(round(y[i,1], 7)!=round(y[i-1,1],7)){
+   y<-rbind(y,NA)
+   i<-i+1
+   k<-k+1
+   for(j in 1:1000){
+ 
+     x<-matrix(dat=c(0, y[i-1,j]*0, 3.880258, 0, 0*0, 0, 0*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
+     y[i,j]<-x%*%t(t(beta.draws[j,]))+8.8
+     #print(y[i,1], quote=F)
+     
+   }
+   
+ }
> 
> 
> y<-rbind(y,NA)
> i<-i+1
> for(j in 1:1000){
+   
+   x<-matrix(dat=c(1, y[i-1,j]*1, 3.880258, 0, 0*1, 0, 0*1, .3926666, .3926666*1, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
+   y[i,j]<-x%*%t(t(beta.draws[j,]))+8.8
+   #print(y[i,1], quote=F)
+   
+ }
> 
> pen<-i+9
> 
> while(i<pen){
+   
+   y<-rbind(y,NA)
+   i<-i+1
+   for(j in 1:1000){
+     
+     x<-matrix(dat=c(1, y[i-1,j]*1, 3.880258, 0, 0*1, 0, 0*1, .3926666, .3926666*1, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
+     
+     y[i,j]<-x%*%t(t(beta.draws[j,]))+8.8
+     #print(y[i,1], quote=F)
+     
+   }
+   
+ }
> 
> 
> while(round(y[i,1], 7)!=round(y[i-1,1],7)){
+   y<-rbind(y,NA)
+   i<-i+1
+   for(j in 1:1000){
+     
+     x<-matrix(dat=c(0, y[i-1,j]*0, 3.880258, 0, 0*0, 0, 0*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
+     y[i,j]<-x%*%t(t(beta.draws[j,]))+8.8
+     #print(y[i,1], quote=F)
+     
+   }
+   
+ }
> 
> end.point<-(k-10)
> t<-0:(i-end.point)
> plot(y[end.point:i,1]~t, type="l", xlab="time")
> abline(v=10, lty=2)
> abline(v=20, lty=2)
> 
> 
> traj.quant<-apply(X=y, MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
> 
> 
> pdf("dyadaidtrajectory.eps")
> library(gplots)
> plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=traj.quant[1,end.point:i], ylab="ln(bilateral aid PC + 1)", xlab="time", main = c("simulated aid trajectory with 90% CI error bars"), ylim=c(5.78, 10.02))
> rect(11,0,20,15,density=NA, col="grey85")
> plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=traj.quant[1,end.point:i], add=T)
> box()
> abline(v=11, lty=2, lwd=2)
> abline(v=20, lty=2, lwd=2)
> legend("topright", legend=c(" median aid prediction", " condemnation period"), pt.cex=c(1,2), pch=c(1,22), pt.bg="gray85")
> dev.off()
pdf 
  2 
> 
> 
> traj.quant[2, (end.point+10):(end.point+20)]
 [1] 6.692112 7.244575 7.606398 7.822297 7.942991 8.038885 8.095074 8.126970
 [9] 8.146867 8.159559 8.167936
> # percentage change in aid
> exp(traj.quant[2, (end.point+20)]) / exp(traj.quant[2, (end.point+10)]) - 1
     50% 
3.374639 
> 
> 
> 
> 
> 
> x<-matrix(dat=c(0, 2*0, 3.880258, 0, 0*0, 0, 0*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, 2, 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
> y<-rbind(1,x%*%t(beta.draws)+6)
> i<-2
> k<-2
> 
> while(round(y[i,1], 7)!=round(y[i-1,1],7)){
+   y<-rbind(y,NA)
+   i<-i+1
+   k<-k+1
+   for(j in 1:1000){
+     
+     x<-matrix(dat=c(0, y[i-1,j]*0, 3.880258, 0, 0*0, 0, 0*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
+     y[i,j]<-x%*%t(t(beta.draws[j,]))+6
+     #print(y[i,1], quote=F)
+     
+   }
+   
+ }
> 
> 
> y<-rbind(y,NA)
> i<-i+1
> for(j in 1:1000){
+   
+   x<-matrix(dat=c(1, y[i-1,j]*1, 3.880258, 0, 0*1, 0, 0*1, .3926666, .3926666*1, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
+   y[i,j]<-x%*%t(t(beta.draws[j,]))+6
+   #print(y[i,1], quote=F)
+   
+ }
> 
> pen<-i+9
> 
> while(i<pen){
+   
+   y<-rbind(y,NA)
+   i<-i+1
+   for(j in 1:1000){
+     
+     x<-matrix(dat=c(1, y[i-1,j]*1, 3.880258, 0, 0*1, 0, 0*1, .3926666, .3926666*1, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
+     y[i,j]<-x%*%t(t(beta.draws[j,]))+6
+     #print(y[i,1], quote=F)
+     
+   }
+   
+ }
> 
> 
> while(round(y[i,1], 7)!=round(y[i-1,1],7)){
+   y<-rbind(y,NA)
+   i<-i+1
+   for(j in 1:1000){
+     
+     x<-matrix(dat=c(0, y[i-1,j]*0, 3.880258, 0, 0*0, 0, 0*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
+     y[i,j]<-x%*%t(t(beta.draws[j,]))+6
+     #print(y[i,1], quote=F)
+     
+   }
+   
+ }
> 
> end.point<-(k-10)
> t<-0:(i-end.point)
> plot(y[end.point:i,1]~t, type="l", xlab="time")
> abline(v=10, lty=2)
> abline(v=20, lty=2)
> 
> 
> traj.quant<-apply(X=y, MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
> 
> 
> library(gplots)
> pdf("dyadaidtrajectory_small.eps")
> plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=pmax(traj.quant[1,end.point:i],0), ylab="ln(bilateral aid PC + 1)", xlab="time", main = c("simulated aid trajectory with 90% CI error bars"))
> rect(11,-1,20,15,density=NA, col="grey85")
> plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=pmax(traj.quant[1,end.point:i],0), add=T)
> box()
> abline(v=11, lty=2, lwd=2)
> abline(v=20, lty=2, lwd=2)
> legend("bottomright", legend=c(" median aid prediction", " condemnation period"), pt.cex=c(1,2), pch=c(1,22), pt.bg="gray85")
> dev.off()
pdf 
  2 
> 
> 
> traj.quant[2, (end.point+10):(end.point+20)]
 [1] 2.1283470 1.5732516 1.2104380 0.9964787 0.8733915 0.7833948 0.7233917
 [8] 0.6851998 0.6652119 0.6507730 0.6410873
> # percentage change in aid
> 1 - exp(traj.quant[2, (end.point+20)]) / exp(traj.quant[2, (end.point+10)])
      50% 
0.7740089 
> 
> 
> 
> 
> 
> 
> 
> 
> x<-matrix(dat=c(0, 2*0, 3.880258, 0, 0*0, 1, 1*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, 2, 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
> y<-rbind(1,x%*%t(beta.draws)+6)
> i<-2
> k<-2
> 
> while(round(y[i,1], 7)!=round(y[i-1,1],7)){
+   y<-rbind(y,NA)
+   i<-i+1
+   k<-k+1
+   for(j in 1:1000){
+     
+     x<-matrix(dat=c(0, y[i-1,j]*0, 3.880258, 0, 0*0, 1, 1*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
+     y[i,j]<-x%*%t(t(beta.draws[j,]))+6
+     #print(y[i,1], quote=F)
+     
+   }
+   
+ }
> 
> 
> y<-rbind(y,NA)
> i<-i+1
> for(j in 1:1000){
+   
+   x<-matrix(dat=c(1, y[i-1,j]*1, 3.880258, 0, 0*1, 1, 1*1, .3926666, .3926666*1, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
+   y[i,j]<-x%*%t(t(beta.draws[j,]))+6
+   #print(y[i,1], quote=F)
+   
+ }
> 
> pen<-i+9
> 
> while(i<pen){
+   
+   y<-rbind(y,NA)
+   i<-i+1
+   for(j in 1:1000){
+     
+     x<-matrix(dat=c(1, y[i-1,j]*1, 3.880258, 0, 0*1, 1, 1*1, .3926666, .3926666*1, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
+     y[i,j]<-x%*%t(t(beta.draws[j,]))+6
+     #print(y[i,1], quote=F)
+     
+   }
+   
+ }
> 
> 
> while(round(y[i,1], 7)!=round(y[i-1,1],7)){
+   y<-rbind(y,NA)
+   i<-i+1
+   for(j in 1:1000){
+     
+     x<-matrix(dat=c(0, y[i-1,j]*0, 3.880258, 0, 0*0, 1, 1*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
+     y[i,j]<-x%*%t(t(beta.draws[j,]))+6
+     #print(y[i,1], quote=F)
+     
+   }
+   
+ }
> 
> end.point<-(k-10)
> t<-0:(i-end.point)
> plot(y[end.point:i,1]~t, type="l", xlab="time")
> abline(v=10, lty=2)
> abline(v=20, lty=2)
> 
> 
> traj.quant<-apply(X=y, MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
> 
> 
> pdf("dyadaidtrajectory_alliance.eps")
> library(gplots)
> plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=pmax(traj.quant[1,end.point:i],0), ylab="ln(bilateral aid PC + 1)", xlab="time", main = c("simulated aid trajectory with 90% CI error bars"))
> rect(11,-1,20,15,density=NA, col="grey85")
> plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=pmax(traj.quant[1,end.point:i],0), add=T)
> box()
> abline(v=11, lty=2, lwd=2)
> abline(v=20, lty=2, lwd=2)
> legend("topright", legend=c(" median aid prediction", " condemnation period"), pt.cex=c(1,2), pch=c(1,22), pt.bg="gray85")
> dev.off()
pdf 
  2 
> 
> 
> traj.quant[2, (end.point+10):(end.point+20)]
 [1] 1.998582 2.230155 2.358357 2.449682 2.509200 2.548075 2.572129 2.590006
 [9] 2.597344 2.602000 2.604130
> # percentage change in aid
> 1 - exp(traj.quant[2, (end.point+20)]) / exp(traj.quant[2, (end.point+10)])
       50% 
-0.8322562 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> ### plot instantaneous MEs, Murdie/Davis Model 2
> 
> 
> rm(list=ls())
> library(MASS)
> require(foreign)
> 
> dat.neilsen <- read.dta(file="ngo-shaming.dta")
> 
> neilsen.beta<-read.delim(file="tobit_beta_ngo.txt", header=T)
> neilsen.beta<-as.matrix(neilsen.beta[2:26])
> names(neilsen.beta)<-NULL
> 
> neilsen.vcv<-read.delim(file="tobit_VCV_ngo.txt")
> neilsen.vcv<-as.matrix(neilsen.vcv[1:25,2:26])
> 
> beta.draws<-mvrnorm(n=1000, mu=neilsen.beta, Sigma=neilsen.vcv)
> 
> 
> # instantaneous effects of NGO Shaming
> par(mar=c(5,5,4,2))
> y.seq<-seq(from=0, to=13, by=0.1)
> me<-matrix(data=NA, ncol=3, nrow=length(y.seq))
> 
> for(i in 1:length(y.seq)){
+   
+   y.test<-y.seq[i]
+   me.draws<-beta.draws[,1]+beta.draws[,2]*y.test
+   me[i,]<-quantile(me.draws, probs=c(0.025, 0.5, 0.975))
+   
+   
+ }
> 
> pdf("ngo_instant.pdf")
> # plot instantaneous marginal effects
> plot(NA, ylim=c(-0.25, 0.25), xlim=c(0,13), xaxt="n", main=c("Instantaneous Marginal Effect","of NGO Shaming on Bilateral Economic Aid"), xlab=expression(ln("Bilateral Economic Aid + 1")[it-1]), ylab=expression(partialdiff*ln("Bilateral Economic Aid + 1")[it]/partialdiff*Shaming[i(t-1)]))
> axis(side=1, at=seq(from=0, to=12, by=2))
> hist.back<-hist(dat.neilsen$l_lneconaidpc, plot=FALSE)
> hist.back$density <- hist.back$density*0.3
> lines(hist.back, border=hsv(0, 0, 0.5, alpha=1), col=hsv(0, 0, 0.5, alpha=0.25), freq=FALSE, alpha=1)
> 
> lines(me[,2]~y.seq)
> lines(me[,1]~y.seq, lty=2)
> lines(me[,3]~y.seq, lty=2)
> abline(h=0, lty=3)
> 
> legend("topright", legend=c("marginal effect", "95% confidence interval"), lty=c(1,2))
> text(x=10, y=-2, "histogram indicates 3*(data density)", cex=0.75)
> dev.off()
pdf 
  2 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> rm(list=ls())
> library(MASS)
> require(foreign)
> 
> dat.neilsen <- read.dta(file="ngo-shaming.dta")
> 
> neilsen.beta<-read.delim(file="tobit_beta_ngo2.txt", header=T)
> neilsen.beta<-as.matrix(neilsen.beta[2:26])
> names(neilsen.beta)<-NULL
> 
> neilsen.vcv<-read.delim(file="tobit_VCV_ngo2.txt")
> neilsen.vcv<-as.matrix(neilsen.vcv[1:25,2:26])
> 
> beta.draws<-mvrnorm(n=1000, mu=neilsen.beta, Sigma=neilsen.vcv)
> 
> 
> # instantaneous effects of NGO Shaming
> par(mar=c(5,5,4,2))
> y.seq<-seq(from=0, to=13, by=0.1)
> me<-matrix(data=NA, ncol=3, nrow=length(y.seq))
> 
> for(i in 1:length(y.seq)){
+   
+   y.test<-y.seq[i]
+   me.draws<-beta.draws[,1]+beta.draws[,2]*y.test
+   me[i,]<-quantile(me.draws, probs=c(0.025, 0.5, 0.975))
+   
+   
+ }
> 
> pdf(file="ngo_dum_instant.pdf")
> # plot instantaneous marginal effects
> plot(NA, ylim=c(-1.5, 1.5), xlim=c(0,13), xaxt="n", main=c(""), xlab=expression(ln("Bilateral Economic Aid + 1")[it-1]), ylab=expression(partialdiff*ln("Bilateral Economic Aid + 1")[it]/partialdiff*Shaming[i(t-1)]))
> mtext(line=0.75, expression(paste("NGO Shaming">=1~"on Bilateral Economic Aid")), cex=1.1)
> mtext(line=2, "Instantaneous Marginal Effect of", cex=1.1)
> axis(side=1, at=seq(from=0, to=12, by=2))
> hist.back<-hist(dat.neilsen$l_lneconaidpc, plot=FALSE)
> hist.back$density <- hist.back$density*2
> lines(hist.back, border=hsv(0, 0, 0.5, alpha=1), col=hsv(0, 0, 0.5, alpha=0.25), freq=FALSE, alpha=1)
> 
> lines(me[,2]~y.seq)
> lines(me[,1]~y.seq, lty=2)
> lines(me[,3]~y.seq, lty=2)
> abline(h=0, lty=3)
> 
> legend("topright", legend=c("marginal effect", "95% confidence interval"), lty=c(1,2))
> text(x=10, y=-2, "histogram indicates 3*(data density)", cex=0.75)
> 
> dev.off()
pdf 
  2 
> 
> 
> proc.time()
   user  system elapsed 
  20.54    0.96   20.56 
